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ABSTRACT 

We present for astrophysical use a multi- dimensional numerical code to solve the 
equations for ideal magnetohydrodynamics (MHD). It is based on an explicit finite 
difference method on an Eulerian grid, called the Total Variation Diminishing (TVD) 
scheme, which is a second-order-accurate extension of the Roe-type upwind scheme. 
Multiple spatial dimensions are treated through a Strang-type operator splitting. The 
constraint of a divergence-free field is enforced exactly by calculating a correction via a 
gauge transformation in each time step. 

Results from two-dimensional shock tube tests show that the code captures cor- 
rectly discontinuities in all three MHD waves families as well as contact discontinu- 
ities. The numerical viscosities and resistivity in the code, which are useful in order to 
understand simulations involving turbulent flows, are estimated through the decay of 
two-dimensional linear waves. Finally, the robustness of the code in two-dimensions is 
demonstrated through calculations of the Kelvin-Helmholtz instability and the Orszag- 
Tang vortex. 

Subject headings: hydromagnetics - magnetohydrodynamics:MHD - methods :numerical 
- shock waves 
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1. INTRODUCTION 

In recent years upwind finite difference schemes based on ideal gas conservation laws have been 
popular for solving compressible hydrodynamic conservation equations, mainly due to their ability 
to sharply capture discontinuities and also to the robustness of the schemes. Several examples 
include Godunov's Scheme (Godunov 1959), MUSCL scheme (Van Leer 1979), Roe Scheme (Roe 
1981), TVD scheme (Harten 1983), PPM scheme (Colella k Woodward 1984), and ENO scheme 
(Harten et al. 1987). A comparative review for some of these schemes can be found in Woodward 
k Colella (1984). Such methods have been very effectively applied to many astrophysical problems 
{e.g., cosmological hydrodynamics with the TVD scheme (Ryu et al. 1993), accretion flow with the 
ENO scheme (Ishii et al. 1993), supernova remnants (Chevalier, Blondin k Emmering 1992) and 
astrophysical convection (Porter k Woodward 1994) with the PPM scheme). 

As important as these hydrodynamical tools have been, there are many applications that 
fundamentally depend on the presence of magnetic fields and an ionized, conducting character of the 
medium under investigation. To the extent that these media can be treated through continuum fluid 
mechanics, this means one must work with magnetohydrodynamic (MHD) equations rather than 
hydrodynamic equations. But, due to the intrinsic complexity of the MHD flows, the development 
of numerical techniques to solve MHD equations has been slower than for hydrodynamics. For 
instance, until now most numerical schemes have been based on methods dependent on artificial 
viscosity to form shocks(e.</., DeVore 1991; Lind, Payne k Meier 1991; Stone k Norman 1992). 
Those schemes have been used successfully in astrophysical applications {e.g., Lind et al. 1989; 
Stone k Norman 1994). However, since past experience with fully conservative, high-order upwind 
hydrodynamic schemes found those to be superior in many applications (Woodward k Colella 1984), 
it is naturally interesting to extend such schemes to solve MHD conservation equations. Several 
investigators who have worked in recent years on the development of high-order, conservative, 
upwind differencing schemes for MHD include Brio k Wu (1988), Zachary k Colella (1992), 
Zachary, Malagoli, k Colella (1994), and Dai k Woodward (1994a,1994b). Brio k Wu applied 
Roe's approach to the MHD equations. Zachary and collaborators used the BCT scheme (Bell, 
Colella, k Tragenstein 1989) to estimate fluxes in MHD conservation equations. Dai k Woodward 
applied the PPM scheme to MHD. One of the special concerns in developing a scheme for MHD 
is the fact that the equations may have solutions at which they are no longer strictly hyperbolic, 
meaning some of the characteristic speeds become locally degenerate. 

In a previous work (Ryu k Jones 1995), we described a one- dimensional conservative numerical 
code to solve the equations for ideal magnetohydrodynamics. It is based on an explicit finite 
difference scheme on an Eulerian grid, called the Total Variation Diminishing (TVD) scheme 
(Harten 1983), which is a second-order- accurate extension of the Roe-type upwind scheme. Tests 
using an extensive set of MHD shock tube problems showed that the one- dimensional code can 
resolve strong shocks within 2-4 cells; weaker shocks, especially slow shocks are somewhat broader. 
Contact and tangential discontinuities also require more cells for capture. 
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This code provides a promising first step to a useful computational tool. However, there 
are only a very limited number of practical applications that can be studied under the planar 
symmetry constraint. Most real flows will depend on structures in at least two or, more likely, in 
three spatial dimensions. Thus, general application of our code requires that it be extended to a 
multi- dimensional form. In this paper, we describe this extension. In the multi- dimensional code, 
multiple spatial dimensions are treated through a Strang-type operator splitting (Strang 1968). 
The constraint of a divergence-free field is enforced exactly by calculating a correction via a gauge 
transformation in each time step (Brackbill & Barnes 1980). 

This scheme should have many interesting application to astrophysical problems. We have 
already applied the one- dimensional version of our code successfully to the first time dependent 
study of diffusive cosmic-ray acceleration in oblique MHD shocks (Frank, Jones & Ryu 1994, 1995a). 
We demonstrated that inclusion of magnetic fields in the fluid dynamics introduces a number of 
effects, including some that are subtle and not apparent except in a time dependent simulation. 
With the new multi- dimensional code described here, we can extend that analysis to include such 
important issues as the evolution of unstable MHD cosmic-ray shocks. Tests we describe below 
illustrate some of the other kinds of problems that can be addressed with the code. We have 
underway a study of the nonlinear evolution of the MHD Kelvin-Helmholtz instability (Frank, 
Jones & Ryu 1995b), for example. 

This paper is organized in the following way. In §2, we describe the extension of the one- 
dimensional code into a multi- dimensional form, including the operator splitting scheme and the 
correction to preserve V B = 0. In §3, we present the results of numerical tests performed with the 
two-dimensional version of the code, including shock tube problems, the decay of linear waves, the 
growth of the Kelvin-Helmholtz instability, and evolution of the Orszag-Tang vortex. Our findings 
are summarized in §4. 



2. THE NUMERICAL SCHEME 



2.1. The MHD Equations 



MHD describes the behavior of the combined system of a conducting fluid and magnetic fields in 
the limit that the displacement current and the separation between ions and electrons are neglected. 
So, the MHD equations represent coupling of the equations of fluid dynamics with the Maxwell's 
equations of electrodynamics. Here, we describe a numerical scheme to calculate the evolution of 
the following ideal MHD equations, where the effects of electrical resistivity, viscosity, and thermal 
conductivity are dropped, 

^ + V-(pv) = 0, (2.1) 
dv 11 

-7T- + v ■ Vv + -Vp- -(V x B) x B = 0, (2.2) 
at p p 
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dp 
di 
dB 



— + v ■ Vp + 7PV • v = 0, 



at 



V x x B) 



0, 



(2.3) 



(2.4) 



(see Jackson (1975) for the derivation of the equations). Here, we have chosen units so that factor 
of 47T does not appear in the equations. An additional explicit constraint V • B = is imposed 
to account for the absence of magnetic monopoles. Although that constraint is included in the 
derivation of equation (2.4), it generally cannot be maintained precisely in differenced forms of 
that equation, even when one works with an exactly conservative scheme as we will outline. 



In Cartesian geometry, the above equations are written in conservative form as 
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with Fy and Fz obtained by properly permuting indices. Here, the total pressure and the total 
energy are given by 
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(2.7) 



(2.8) 



With the state vector, q, and the flux functions, Fx(q), Fy(q), and Fz(q), the Jacobian matrices, 
A x (q) = dFx/dq, Ay(q) = dFy/dq, and A z (q) = dF z /dq, are formed. The system of 
equations is called hyperbolic if all the eigenvalues of the Jacobian matrices are real and distinct 
and the corresponding set of right eigenvectors is complete (Jeffrey & Taniuti 1964). The MHD 
equations form a non- strictly hyperbolic system, meaning that some eigenvalues may coincide at 
some points. 



2.2. The One-Dimensional Code 



The procedure to build a one- dimensional MHD code based on the TVD scheme was described 
in detail in our previous paper (Ryu & Jones 1995). Here, we briefly summarize it and mention 
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optimizations for the multi- dimensional extension. To start, we consider a plane-symmetric, one- 
dimensional flow exhibiting variation along the x-direction. The first step to build the code is to 
find the eigenvalues and the right and left eigenvectors of the Jacobian matrix, Ax(q)- The seven 
eigenvalues a\ , • • • , in non-increasing order are 



ai ,7 = v x ± c 



a 2 ,e = v x ± c a , 



a 3 , 5 = v x ± c s , 



a 4 = v X: 



(2.9) 



where cy, Ca, c$ are the fast, Alfven, and slow characteristic speeds. The quantities ai,---,ar 
represent the seven speeds with which information is propagated locally by three MHD wave families 
and an entropy mode. The three characteristic speeds are expressed as 



7,* 




(2.10) 



and 



(2.11) 



with the sound speed, 




(2.12) 



The corresponding eigenvectors are given, for example, in Jeffrey & Taniuti (1964). However, 
near a point where either Bx = or By = Bz = 0, the eigenvectors are not well defined, with 
elements becoming singular. By renormalizing the eigenvectors, the singularities can be removed. 
The renormalized eigenvectors are given in Brio & Wu (1988) and Ryu & Jones (1995). 

In a code based on the TVD scheme, the physical quantities are referred to the cell centers while 
the fluxes are computed on the cell boundaries. Implementation of Roe's linearization technique 
would result in a particular averaged form of the physical quantities on the cell boundaries (Roe 
1981). However, as pointed out by Brio & Wu (1988), it is not possible to derive this particular 
analytic form of the averaged quantities in MHD for general cases with an adiabatic index 7^2. 



Instead, we modify Roe's scheme and use p 
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on the cell boundaries with the arithmetic averages at i and i + 1. Then, other quantities like 
momentum, gas pressure, total energy, etc are calculated by combining those quantities. Our tests 
for cases with 7 = 2 indicated that the above simple averaging would do just as well when compared 
to the full implementation of Roe's linearization technique, as already suggested for the original 
hydrodynamic code by Harten (1983). 

The state vector q n at the time step n is updated by calculating the modified flux fx at the 
cell boundaries as follows: 
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2.15) 
2.16) 

2.17) 

2.18) 
2.19) 

2.20) 
2.21) 



Here, a^, Rfc, are the seven eigenvalues in (2.9) and the corresponding right and left 
eigenvectors. The time step At n is restricted by the usual Courant condition for the stability, 



At n = C COUT Ax/Msix(\v 



n 

X *Z-\--^ 



+ c 



n 



with C cour < 1. Typically we use C c 



0.6, although 



values up to 0.9 seems to be sufficient for most calculations. Note that the e and C cour values in 
(2.20) as optimized for the multi- dimensional code are different from those recommended for the 
one- dimensional code (Ryu & Jones 1995). Our tests showed that procedures to steepen contact 
discontinuities and rotational discontinuities similar to those suggested in the original TVD paper 
by Harten (1983) do not work very well in this context. They produce spurious numerical oscil- 
lations that can significantly degrade the solution. So, we do not include contact steepener and 
rotational steepener routines in our multi- dimensional MHD code. 



2.3. The Multi-Dimensional Extension 



The multi- dimensional extension has been done through a Strang-type directional splitting 
(Strang 1968). In each time step, multi- dimensional derivatives are split into a set of one- 
dimensional derivatives, with variations in other directions ignored temporarily. Then, each row 
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and column in the grid is treated as if it were a one- dimensional problem. Updating the flow 
quantities along each row is done using the one- dimensional code described in the previous section. 
The parallel (to the direction of the row) component of magnetic field is kept constant and only 
the perpendicular component is updated. One complete time step updating the full state vector 
q n to q n+1 in each cell is composed of updating it along two or three directions, as appropriate. 
For instance, in three-dimensional Cartesian geometry, the state vector is updated along x, y, and 
^-directions, so 

q n+1 = L z L y L x q n . (2.22) 

In order to maintain second-order accuracy, the order of directional passes is permuted by 
the Strang-type prescription: L z L y L x , L x L y L z , L x L z L y , L y L z L x , L y L x L z , and then 
L z L x L y , for example. The time step is restricted to satisfy the Courant condition along each 
row in three directions. It is calculated at the start of the above permuting sequence and used 
through one complete sequence. 



2.4. Ensuring V B = 



The condition V • B = is, of course, a necessary initial constraint in multi- dimensional 
MHD flows and should be preserved during their evolution. While the differential MHD equations 
formally ensure V ■ B = 0, the numerical errors due to discretization and operator splitting can lead 
to non-zero divergence over time. Physically, this is due to the fact that even though schemes such 
as this one are exactly conservative of the fluxes in equation (2.6) on the cell boundaries crossing 
each directional pass, nothing maintains a constant magnetic flux in the Gauss's Law sense across 
the entire surface of each cell. That is, nothing forces conservation of magnetic charge over a finite 
cell during the entire time step. This error usually grows exponentially during the computations, 
causing an artificial force parallel to the magnetic field and destroying the correct dynamics of the 
flows (see e.g., Schmidt- Voigt 1987). 

Zachary, Malagoli, & Colella (1994) described that, by using the modified form of the MHD 
equations originally suggested in Brackbill & Barnes (1980), the error in V • B could be kept 
small enough in their code that their tests showed no discernible differences with and without the 
correction for non-zero V • B. So they abandoned the correction in their code. However, since the 
modified form is not suitable for our code and it does not remove V • B completely anyway, we 
enforce V • B = explicitly by a simple transformation. First we solve for the potential, </>, defined 
by the Poisson equation 

V 2 </» + V • B = 0. (2.23) 

Here, B is the updated magnetic field obtained by the procedure already outlined. Then we 
compute the corrected magnetic field, defined as B c = B + V</>, for which V • B c = 0, as required. 
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The difference form of the above Poisson equation in two-dimensional Cartesian geometry is 



4 (Ax) 2 4 (Ay) 2 

= _ f B x,i+i,j ~ B x,i-i,j + B y,i,j+i ~ ^y,zj-A ^ ^ 

To solve the above difference equation, we employ a suitably selected rapid technique. On a 

Cartesian grid, for example, we can take advantage of the fact from the Fourier transform derivative 

theorem that FFTs (fast Fourier transforms) can provide an exact solution to </> ■ ■ from equation 

'v. 

(2.23) if the field structure is periodic on some space (e.g., Bracewell 1986). Even if the magnetic 
field is not periodic on the computational box of interest, it is often possible to create such a 
structure for this transformation by doubling the box size for this step only and choosing an 
appropriate symmetry within the extended space. According to our tests, the cost to enforce 
V • B = explicitly is reasonable, taking ^5% of the total CPU time in cases with periodic 
boundaries on the computational space. With others, like reflecting or continuous boundaries, that 
require an extended space the cost to enforce V • B = is ^30% of the total CPU time. 

The present version of the multi- dimensional MHD TVD code runs at ~ 400 MFlops on a single 
Cray C90 processor, although it is not specifically optimized for the machine. This corresponds to 
an update rate of ~ 2 X 10 5 cells/sec for the one- dimensional version and ~ 10 5 cells/sec for the 
two-dimensional version. 



3. TESTS 



In this section we briefly describe several tests we have carried out to determine the charac- 
teristics of the code and its suitability for practical application. We aim to examine the ability of 
the code to handle all three MHD wave families as well as its performance in computing complex 
flows. In addition we will estimate the effects of numerical viscosity and electrical resistivity. For 
each of these tests we assume an adiabatic index, 7 = 5/3 and a Courant constant, C'cour = 0.6. 



3.1. Shock Tube tests 



Extensive tests of the code have been done with MHD shock tube problems, using a two- 
dimensional Cartesian grid and various orientations of the structures with respect to the grid. 
Here we describe tests performed in a two-dimensional box with x = [0, 1] and y = [0, 1], where 
structures like discontinuities and waves propagate along the diagonal line joining x = y = and 
x = y = 1. We present two tests. One includes only two (x and y) components of magnetic field 
and velocity, so that they are confined in the computational plane. The second includes all three 
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components. The results from the numerical calculations are marked with dots in Figures 1 and 
2. They can be compared with analytic solutions from the nonlinear Riemann solver described in 
Ryu and Jones (1995), plotted here with lines. Structures are measured along the diagonal line 
joining x = y = and x = y = 1. The plotted quantities are density, gas pressure, total energy, v\\ 
(velocity parallel to the diagonal line; i.e., parallel to the wave normal), Vj_ (velocity perpendicular 
to the diagonal line but still in the computational plane), Vz (velocity in the direction out of plane), 
and the analogous magnetic field components, B\\, B_\_, and Bz- 

The first test in Fig. 1 has been done with two-dimensional magnetic field and velocity vectors 
in the x — y plane. In this case the sign of the perpendicular magnetic field (Bj_) remains unchanged 
across the structures. The initial left state is (p, v\\, Vj_, Vz, B±, Bz, E) = (1, 10, 0, 0, 5/a/47T, 
0, 20) and the initial right state is (1, -10, 0, 0, 5/\/4tt, 0, 1), with B\\ = 5/V4tt. Fig. la in Ryu 
& Jones (1995) illustrates a solution to this problem with the one- dimensional MHD TVD code. 
The calculation has been done using 256 X 256 cells, and plots correspond to time t = 0.08\/2- The 
structure is bounded by a left and right facing fast shock pair. There are also a left facing slow 
rarefaction, a right facing slow shock and a contact discontinuity. Three anomalous points visible 
in the Vj_ plot are due to the error induced in the subtraction of two big numbers to get a small 
number. Note that v± has been calculated with v± = { — Vx +Vy)/\/2 and Vx and Vy are, at least, 
an order of magnitude larger than Vj_ in preshock regions. As expected from the two-dimensional 
nature of the magnetic field and velocity structure, there is no rotational discontinuity. In this test 
the fast shocks are strong, with a large parallel velocity jump, [v\\]. The density, parallel velocity 
and pressure jumps in the slow shock are weak, so that this feature is mainly apparent through the 
jumps in the tangential velocity, [vj_] and tangential magnetic field, [Bj_]. The capture of shocks 
and the contact discontinuity here are virtually the same as with the one- dimensional code. 

The second test in Fig. 2 involves all three components of both magnetic field and velocity, 
with the magnetic field plane rotated across the initial discontinuity (such calculations are often 
described as 2 + 1/2 dimensional). This test is also presented as Fig. 2a in Ryu & Jones (1995). 
The initial left state is (p, v\\, v ± , v z , B ± , B z , E) = (1.08, 1.2, 0.01, 0.5, 3.6/v^vr, 2/v^vr, 
0.95) and the initial right state is (1, 0, 0, 0, 4/V47r, 2/V47r, 1), with B\\ = 2/\/4tt. Again the 
calculation has been done using 256x256 cells, and plots correspond to time t = 0.2\/2. Fast shocks, 
rotational discontinuities, and slow shocks propagate from each side of the contact discontinuity. 
Here, the rotation of the magnetic field across the initial discontinuity generates two rotational 
discontinuities. Again, the structures are captured by the two-dimensional code very similarly to 
what we found with the one- dimensional code. 

We also carried out many of the other shock tube tests described in Ryu & Jones with this code, 
including rarefactions of both fast and slow wave families. Generally the solutions are comparable 
to those seen there. In summary, then, results from multi- dimensional shock tube tests show that 
the code captures correctly discontinuities and rarefaction waves in all three MHD wave families, 
as well as that carried with the entropy mode. 
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3.2. Decay of Linear Waves 



Our conservation equations (2.5-2.6) refer to ideal MHD without the effects of viscosity, 
electrical resistivity, and thermal conductivity. However, in numerical calculations on a discrete 
grid, diffusion of magnetic field and momentum across cell boundaries is unavoidable and introduces 
effective numerical electrical resistivity and viscosity. Energy diffusion also occurs and produces 
numerical thermal conductivity. But, the effects of the numerical thermal conductivity are mostly 
small compared to those of the numerical viscosity and electrical resistivity. 

In order to estimate the properties of the numerical resistivity and viscosity in the multi- 
dimensional MHD code, we have followed the decay of two-dimensional Alfven, fast, and slow 
waves in numerical calculations and compared their decay rates to the predicted rates in a viscous, 
resistive fluid. The MHD equations for a viscous, resistive fluid can be written as 

^ + V-(pv) = 0, (3.1) 
dv 1 1 1 

— + v.W + -Vp--(VxB)xB = -d k a ik , (3.2) 



dp 
dt 



v -Vp + ^pV-v = (7- l)(7 lk d k v v (3.3) 
dB 



dt 



V x (v x B) = r]V 2 B. (3.4) 



In these equations, 



a ik=l 1 \ d k v i + d i v k~ ? S ik V ■ v )+ ( S ik V ■ v ( 3 - 5 ) 



2 

— < 

3 

is the viscosity tensor, where p and (* are the dynamic shear and bulk viscosity, respectively, i] is 
the electrical resistivity, and & k = d/dx k . 

In the test of the decay of Alfven (shear) waves, we have used a standing wave formed along 
the grid diagonal with initial conditions 



8v z = v amp c a sm [k x x + k y y^j 



Sp = Sp = Sv x = Sv y = SB X = SB y = SB Z = 0, (3.6) 

in a stationary background with po = 1, Po = 1 5 an( i B = Box with Bo = 1. This gives 
characteristic speeds, a = 1.291 and Cg i = 0.7071. The calculations have been done in a square 
periodic box with size L = 1 using 8 X 8, 16 X 16, 32 X 32, 64 X 64, 128 X 128, and 256 X 256 
cells. The x and y-wavenumbers have been set to kx = ky = 2ir/L, so the total wavenumber 
k = ^Jk x + ky = \Z2(27r/_L), and the initial peak amplitude has been set to u amp = 0.1. The 
predicted complex frequency of the wave is 



co = l - ( JL + n) k 2 ± cakjl - -^j- (— -r]) 2 k 2 (3.7) 



2 \po J V 4c a V/°o 
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so the decay rate is 
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Ra = c aL 



The effect of the shear viscosity and the electrical resistivity on Alfvenic turbulence of scale, L, is 
measured through a Reynolds number defined by 

When the Reynolds number corresponding to L is large the fluid should be able to maintain on 
that scale a turbulent flow whose properties are not controlled by the properties of the viscosity 
and resistivity. 

In a similar fashion we have set up tests of the decay of compressive waves. Each involves 
standing waves formed on the grid diagonal. For the fast wave we used initial conditions 
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The same background configuration and computational box have been used as in the case of 
Alfven waves. For comparison, the characteristic wave speeds are c j = 1-518 and c$ = 0.6013. 
In compressive waves with a large wave amplitude, nonlinear effects become important (e.g., 
steepening), so a smaller initial peak amplitude, v amp = 10~ 4 , has been used. The derivation 
of the predicted frequency and decay rate of the fast and slow waves in a viscous, resistive fluid is 
complicated. The analytic expressions can be obtained only up to the first order of p, (*, and rj, in 
the limit 
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a condition usually satisfied. The resulting decay rate for fast waves is 



r 



k 2 



f 



7 



4 C-- 

J V3 po 



Po 



n 



Blp_ 

Po Po 



1 p 
3 po 



Po 



p_ 

Po 



(3.12) 



f3.13) 



while the decay rate for slow waves is 
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As in the case of Alfven waves, we define Reynolds numbers 
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R 
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(3.15) 



for the decay of fast waves 



Rs 



8tt 2 c s 
LT S 



(3.16) 



and for the decay of slow waves. 

Fig. 3 shows an example from a 32 X 32 cell test of the decay of a linear wave, an Alfven wave 
in this case. In most calculations, including this one, the decay pattern fits well to an exponential 
form, giving well-defined decay rates. 

In Fig. 4, the normalized decay rates and Reynolds numbers in the tests using Alfven, fast, 
and slow waves are plotted as a function of the number of cells spanning the length L. Good fits of 
the form T cx and R oc n 2 ells are possible, confirming that our code has second-order accuracy. 

Ryu & Goodman (1994) carried out analogous tests for the hydrodynamic TVD code by examining 
decay of two-dimensional sound waves and plane shear flows (Ryu & Goodman 1994). Those tests 
also produce a numerical Reynolds number that scales as R oc n^ells- ^ n the case with a relatively 
strong field examined here (the plasma /3 = po/(B /2) = 2), the numerical Reynolds number that 
applies to a particular number of cells is, on average, roughly 20 -30 times smaller than that in the 
hydrodynamic TVD code. 

However, the amount of numerical dissipation also depends on background configuration, 
especially on background field strength in MHD codes. So we have followed the decay of a fast wave 
in a background with different magnetic field strength (Bo)- The fast wave becomes a sound wave, 
if B = 0. We have used B = 10, 3, 1, 0.3, 0.1, 0.03, 0.01 0.003, 0.001, and (so (3 = 2 x 10" 2 , 
• • •, 2 X 10 6 , and oo) with 128 X 128 cells. For Bo = 0, the hydrodynamic TVD code has been 
used (Harten 1983; Ryu and Goodman 1994). Except Bo, the background configuration is same 
as that used in the test in Fig. 4. Fig. 5 shows the resulting normalized decay rate and Reynolds 
number as functions of /3. Here, the normalized decay rate is biggest around /3 = 1, which is the 
case of the test in Fig. 4. The reason why the normalized decay rate is decreased with decreasing 
/3 for /3^1 is because the fast wave speed increases faster than the decrease in the decay rate. The 
decay rate in the MHD TVD code converges as /3 goes oo. But still it is significantly larger than 
that in the hydrodynamic TVD code. This is because the MHD TVD code does not reduce exactly 
to the hydrodynamic TVD code even B = 0. There are 7 characteristic fields in MHD, and 5 
characteristic fields for hydrodynamics are recovered by combining them properly, as was described 
in Ryu and Jones (1995). This procedure puts in some additional numerical diffusion. 

The above tests mean that to reach an inertial range simulations of MHD turbulence, especially 
with strong field (/3 ~ 1), would need to resolve fine structures with more cells than those of 
hydrodynamic turbulence. That could be an important aspect of such calculations. We suspect that 
this will be a common property of many numerical schemes, but are not aware of any previously 
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published tests that address the issue. Computations primarily concerned with MHD shocks, 
however, can be conducted with similar resolutions as those for hydrodynamic shocks. Numerical 
dissipation measured in this subsection influences smooth flows only, and shocks, especially strong 
shocks, are resolved in the MHD TVD code as sharply as in the hydrodynamic TVD code. 

3.3. The Kelvin-Helmholtz Instability 

To evaluate its potential for practical problems, it is essential to test a numerical method's 
ability to simulate complex, fully two-dimensional flows. We present two sample test problems 
that include some challenging properties, such as strong shear and obliquely intersecting shocks. 
The first test involves an unstable, compressible shear layer; i.e., the compressible MHD Kelvin- 
Helmholtz instability. To make this test as clean and decisive as possible, we have chosen as initial 
conditions a smoothly sheared flow with linear, eigenmode perturbations that have precisely known 
solutions (Miura k Pritchett 1982). 

Initial equilibrium configurations exist in which the magnetic field lines lie in planes of constant 
velocity. For perturbations of this with associated wave vectors parallel to the equilibrium velocity 
field, the dynamical role of the magnetic field depends on the angle between the field and the wave 
vector. If the field is perpendicular to the wave vector (that is, out of the computational plane in 
a two-dimensional simulation) the field can only be compressed in perturbations, not stretched, so 
the perturbation behaves qualitatively like a gasdynamic flow. A much more interesting case for 
our purposes begins with the magnetic field parallel to the flow, since perturbations then bend and 
stretch field lines. We have actually tested both cases, but present here only the latter case. 

Following Miura & Pritchett (1982), we can describe a perturbed equilibrium state as: 

p(x,y,t) = po + 8p(x,y,t), 
v x (x,y,t) = V xo (y) + 8v x (x,y,t), 

vy(x,y,t) = 8v y (x,y,t), 
B x (x,y,t) = B xo + 8B x (x,y,t), (3.17) 
By(x,y,t) = 8B y (x,y,t), 
p(x,y,t) = po + 8p(x,y,t), 
v z (x,y,t) = B z (x,y,t) = 0, 
where the equilibrium velocity profile V X o(y) is given by 

Vn , ,y — 0.5-L, 
V xo (y) = tanh(^ ), (3.18) 

and the equilibrium gas density and pressure along with the magnetic field are uniform. The length, 
h, measures the thickness of the shear layer. 
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The flow is defined in a square computational box with x = [0,L], y = [0,L]. The box is 
assumed to be periodic in the x-direction, and has reflecting y-boundaries. Each of the perturbed 
quantities takes a form 

Sf(x, y, t) = f{y) exp(ik x x + iut), (3.19) 

where kx = 2ir / L and f(y) is a complex function found by numerically integrating the linearized 
versions of equations (2.1-2.4) (see Miura & Pritchett (1982)). In practice it is most convenient to 
define these quantities in terms of the perturbed total pressure, 8p*. To obtain that quantity we 
integrated from one of the y-boundaries to the midplane and assumed antisymmetric and symmetric 
continuations through the midplane for the real and imaginary parts, respectively. Further details 
of the method, including boundary conditions can be obtained from the Miura & Pritchett paper. 
From symmetry the real part of the frequency to is zero in the computational reference frame. The 
imaginary part, or growth rate, T, can be obtained by iteration on the solution. For our purposes 
it was sufficient to use values for T obtained in this fashion and published in graphical form by 
Miura k Pritchett. 

Our test assumes po = 1.0; po = 0.6 (sound speed, a = 1); Bx = 0.4 (Alfven speed, 
c a = 0.4); Vo = 1.0. Thus, the sonic Mach number of the flow is Ms = 1 and the Alfvenic 
Mach number is Mj^ = 2.5. We chose h = _L/(87r), which places the y boundaries far enough 
from the strongly sheared flow to produce minimal boundary effects during the linear phase of 
the instability. The linear growth rate for this mode is T = O.68V0/L. For the test shown we 
began with a perturbation such that \Sp*(y = L)\ = 0.001p£). That carries peak perturbations 
\6p*\ = 0.058p£> and \6v y \ = 0.07V o - 

Figure 5 illustrates the early growth of the perturbation as measured by the root mean squared 
transverse velocity, (^y) 1 / 2 , normalized to the initial value. Results determined from simulations 
with three resolutions, 128 X 128, 256 X 256, and 512 X 512 cells are shown along with the predicted, 
exponential growth curve. After an initial transient, the simulated instability grows at the expected 
rate for t ~ l/T, before saturation becomes significant. The transverse velocity saturates at larger 
values in the higher resolution calculations presumably because the effects of numerical dissipation 
are reduced. 

In Figure 6 we present images of the 512 X 512 simulation at t = 1.4/T. Dynamically this time 
is close to the earlier time shown in Fig. 4 of Miura (1984) for a simulation of the same problem. 
Flow properties are apparently very similar although ours is of higher resolution and shows much 
finer detail. The quantities shown in Fig. 6, here, are (from left to right, top to bottom): p; pm 
(magnetic pressure); p (gas pressure); magnetic field lines (contours of vector potential Az). The 
simulation figure has captured the development of an intense magnetic field region left of center on 
the grid and shows clearly how gas is excluded from that structure. Not obvious in this view, but 
apparent in displays of the flow velocity is a nascent vortex just to the right of center. 



3.4. The Orszag-Tang Vortex 
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In the final test, we have followed the formation of the compressible Orszag-Tang vortex. 
The problem was originally studied by Orszag & Tang (1979) in the context of incompressible 
MHD turbulence, and was later extended by Dahlburg & Picone (1989) and Picone & Dahlburg 
(1991) to compressible MHD turbulence. Zachary, Malagoli, & Colella (1994) applied this problem 
to test their code. We have used this problem to demonstrate the robustness of our code in 
handling problems involving multiple shocks and their collisions and also to compare our code, 
at least qualitatively, with the code by Zachary, Malagoli, & Colella (1994) in a complex, albeit 
two-dimensional, MHD flow. 

The test has been set up with the initial velocity and magnetic field given by 

v = v [— sin(27ry)a; + sin(27rx)y] , (3.20) 

B = B [- sin(27ry)a; + sm{ATtx)y] , (3.21) 

with Vo = 1 and Bo = l/yf&K. Both the velocity and magnetic field contain similar x-points, but 
they have different modal structures. The background density and pressure have been assumed to 
be uniform with values fixed by 

M 2 ee V ° = 1, (3.22) 
(iPo/po) 

= Po - — (3 23) 

with 7 = 5/3. The calculation has been done in a periodic box with x = [0, 1] and y = [0, 1] using 
256 X 256 cells. 

Fig. 7 shows the gray scale images of gas pressure, magnetic pressure, as well as compression, 

V • v, and vorticity, (V X v) z at time t = 0.48. The overall flow structure indicates that the gas 
pressure and the magnetic pressure are anti-correlated resulting in the smoother total pressure, 
except in the postshock flows where both the gas and magnetic pressures increase sharply. The 

V • v plot traces the loci of many shocks, which are interacting with each other. Usually shocks, 
especially strong shocks, are resolved sharply within 2 — 3 cells. The (V X v) z plot demonstrates 
the existence of many fine sheared structures that are not visible in other plots. Even though we 
have neither used exactly the same initial conditions nor made plots at exactly the same epoch, the 
images show that the overall shape and dynamics match closely with those in Dahlburg & Picone 
(1989) and Zachary, Malagoli, k Colella (1994). 



4. Summary 



In this paper, we have described a multi- dimensional numerical code to solve the ideal MHD 
equations. It is based on an explicit and conservative finite difference scheme on an Eulerian grid, 
called the TVD scheme (Harten 1983). The TVD scheme is a second-order- accurate extension 
of the Roe's upwind scheme (Roe 1981). We previously developed and tested a one- dimensional 
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version of the code (Ryu & Jones 1995) and carried out some successful applications, as described 
in §1. The multi- dimensional extension has been done through a Strang-type directional splitting 
(Strang 1968). The constraint V • B = has been enforced exactly by calculating a correction via 
a gauge transformation at each time step. 

We have carried out a number of tests of the method using a two-dimensional Cartesian grid to 
demonstrate its robustness. From shock tube tests we find that the method successfully captures all 
types of MHD discontinuities. It captures strong fast mode shocks generally within 2-3 numerical 
cells, while weak, slow mode shocks may spread over more cells. Other discontinuities, such as 
rotational discontinuity, contact discontinuity, and tangential discontinuity, generally require more 
cells for containment too. We tested the diffusive and dissipative properties of the code by examining 
the decay of three types of linear MHD waves. From this we conclude that to bring turbulent flows 
within an inertial, non-dissipative range requires that important structures exceed ~ 50 cells. To 
examine the ability of the code to correctly follow complex multi- dimensional MHD flows we carried 
out simulations of the evolution of an unstable shear flow and also followed the evolution of the 
Orszag-Tang vortex. In both cases our results agree with previously published behaviors. 

The method is reasonably fast. Our current version of the code runs at about 400 MFlops on 
a Cray C90 processor, so that about 10 5 two-dimensional zones can be updated per second. 

These properties seem to make the code an attractive method for solving practical problems 
involving MHD flows in astrophysical contexts. It is reasonably accurate and robust with small 
enough numerical diffusivity that meaningful solutions to a variety of two-dimensional problems 
ought to be achievable with it using only modest computational resources by today's standards. 
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Research Foundation 1993 and the Basic Science Research Institute Program, Korean Ministry of 
Education 1994, Project No. BSRI-94-5408. TWJ and AF were supported in part by the NSF 
through grants AST91-00486 and AST93-18959, by NASA through grant NAGW-2548 and by the 
University of Minnesota Supercomputer Institute. 



- 17 - 



REFERENCES 

Bell, J. B., Colella, P., k Tragenstein, J. A., 1989, J. of Comput. Phys., 82, 362. 

Bracewell, R. N., 1986, The Fourier Transform and Its Application, (New York: McGraw-Hill). 

Brackbill J. U., k Barnes, D. C, 1980, J. Comp. Phys., 35, 426. 

Brio, M., k Wu, C. C, 1988, J. Comp. Phys., 75, 400. 

Chevalier, R. A., Blondin, J. M. k Emmering, R. T., 1992, ApJ, 392, 118. 

Colella, P., k Woodward, P. R., 1984, J. Comp. Phys., 54, 174. 

Dahlburg, R. B., k Picone, J. M., 1989, Phys. of Fluids B, 1, 2153. 

Dai W., k Woodward P. R., 1994a, J. of Comput. Phys., Ill, 354. 

Dai W., k Woodward P. R., 1994b, preprint. 

DeVore, C. R., 1991, J. of Comput. Phys., 92, 142. 

Frank, A., Jones, T. W., k Ryu, D., 1994, ApJS, 90, 975. 

Frank, A., Jones, T. W., k Ryu, D., 1995a, ApJ, 441, 629. 

Frank, A., Jones, T. W., k Ryu, D., 1995b, in preparation. 

Godunov, S. K., 1959, Math. Sb., 47, 271. 

Harten, A., 1983, J. Comp. Phys., 49, 357. 

Harten, A., Engquist, B., Osher, S., Chakravarthy, S. R., 1987, J. Comp. Phys., 71, 231. 
Ishii, T., Matsuda, T., Shima, E., Livio, M., Anzer, U., k Borner, G., 1993, ApJ, 404, 706. 
Jackson, 1975, Classical Electrodynamics, (New York: John Wiley k Sons, Inc.). 
Jeffrey A., k Taniuti, T., 1964, Nonlinear Waves Propagation, (New York: Academic Press). 
Lind, K. R., Payne, D. G., k Meier, D. L., 1991, preprint. 

Lind, K. R., Payne, D. G., Meier, D. L., k Blandford, R. D., 1989, ApJ, 344, 89. 
Orszag, S. A., k Tang, C. M., 1979, J. Fluid Mech., 90, 129. 
Miura, A. k Pritchett, P. L., 1982, JGR, 87, 7431. 
Miura, A., 1984, JGR, 89, 801. 

Picone, J. M., k Dahlburg, R. B., 1991, Phys. of Fluids B, 3, 29. 

Porter, D. H., k Woodward, P. R., 1994, ApJS, 93, 309. 

Roe, P. L., 1981, J. Comp. Phys., 43, 357. 

Ryu, D., k Goodman, J., 1994, ApJ, 422, 269. 

Ryu, D., k Jones, T. W., 1995, ApJ, 442, 228. 

Ryu, D., Ostriker, J. P., Rang, H., k Cen R., 1993, ApJ, 414, 1. 



- 18 - 



Schmidt- Voigt, M., 1987, in Interstellar Magnetic Fields, ed. Beck, R., k Grave, R., (Berlin: 
Springer), p251. 

Strang, G., 1968, Siam. J. of Num. Anal., 5, 505. 

Stone, J. M., k Norman, M. L., 1992, ApJS, 80, 791. 

Stone, J. M., k Norman, M. L., 1994, ApJ, 433, 746. 

Van Leer, B., 1979, J. Comput. Phys., 32, 101. 

Woodward, P. R., k Colella, P., 1984, J. Comp. Phys., 54, 115. 

Zachary A. L., k Colella P., 1992, J. of Comput. Phys., 99, 341. 

Zachary A. L., Malagoli, A., k Colella P., 1994, Siam. J. Sci. Comp., 15, 263. 



This preprint was prepared with the A AS IATgX macros v3.0. 



- 19 - 



FIGURE CAPTIONS 



Fig. 1 The 2- dimensional MHD shock tube test. Structures propagate diagonally along the line 
(0,0) to (1,1) in the x-y plane. The initial left state is (p, U||, v±, Vz, B±, B z , E) = 
(1, 10, 0, 0, 5/v^Lvr, 0, 20) and the initial right state is (1, -10, 0, 0, 5/\/4tt, 0, 1), with 
= 5/a/Ttt and 7 = 5/3 (same test as Fig. la in Ryu & Jones 1995). Dots are from a 
numerical calculation with the MHD TVD code described in §2 using 256 X 256 cells and a 
Courant constant of 0.6. The structure is shown at time t = 0.08\/2 along the grid diagonal. 
They are plotted along the diagonal line (0, 0) to (1, 1). Lines are from the nonlinear Riemann 
solver described in Ryu & Jones (1995). Plots show from left to right (1) fast shock, (2) slow 
rarefaction, (3) contact discontinuity, (4) slow shock, and (5) fast shock. 

Fig. 2 The 2 + 1/2- dimensional MHD shock tube test. Structures propagate diagonally along the 
line (0,0) to (1, 1) in the x-y plane. The initial left state is (p, v\\, Vj_, Vz, Bz, E) = 
(1.08, 1.2, 0.01, 0.5, 3.6/V4vr, 2/v^vr, 0.95) and the initial right state is (1, 0, 0, 0, 4/v^vr, 
2/\/4tv, 1), with B\\ = 2/V4tt and 7 = 5/3 (same test as Fig. 2a in Ryu & Jones 1995). Dots 
are from a numerical calculation with the MHD TVD code described in §2 using 256 X 256 cells 
and a Courant constant of 0.6. Structure is shown at time t = 0.2\/2, along the diagonal line 
(0, 0) to (1, 1). Lines are from the nonlinear Riemann solver described in Ryu & Jones (1995). 
Plots show from left to right: (1) fast shock, (2) rotational discontinuity, (3) slow shock, (4) 
contact discontinuity, (5) slow shock, (6) rotational discontinuity, and (7) fast shock. 

Fig. 3 The time evolution of the spatially root-mean-squared ^-magnetic field, < SB\ > x / 2 , and 
^-velocity, < 8v 2 z > x / 2 , in the test of the decay of a two-dimensional Alfven wave. An initial 
Alfven wave has been set up in a stationary background and the decay has been followed in 
a numerical calculation using 32 X 32 cells. See text for details. 

Fig. 4 The normalized decay rates, T<jL/ca, T jL/cj, and TgL/cs, and the Reynolds numbers, 
R (see text for definition), versus resolution (the number of cells spanned by the box size L) 
in the tests of the decay of two-dimensional Alfven, fast, and slow waves. Standing waves have 
been set up in a stationary background and the decay rates have been measured during their 
oscillation. The calculations have been done using 8 X 8, 16 X 16, 32 X 32, 64 X 64, 128 X 128, 
and 256 X 256 cells and a Courant constant of 0.6. See text for details. For comparison, 
dotted lines with (YL/cg) oc and R cx n 2 ell are drawn. 

Fig. 5 The normalized decay rate, LL/cy, and the Reynolds number, R, versus the plasma /3 in 
the test of the decay of a two-dimensional fast wave. A standing wave has been set up in 
a stationary background and the decay rate has been measured during its oscillation. The 
calculations have been done with the MHD TVD code for B = 10, 3, 1, 0.3, 0.1, 0.03, 0.01 
0.003, and 0.001. 128 X 128 cells has been used with a Courant constant of 0.6. For Bo = 0, 
the calculation has been done with the hydrodynamic TVD code and its result is indicated 
with dotted lines. See text for details. 
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Fig. 6 The evolution of the spatially root-mean-squared transverse velocity < Vy > x / 2 in the 
MHD Kelvin Helmholtz instability. The results of simulations at three different resolutions are 
shown: (dash-dot = 128 X 128); (dashed = 256 X 256); (solid = 512 X 512). Initial conditions, 
which involve eigenmode perturbations on a smoothly sheared flow of uniform density and 
pressure are described in the text. Time is expressed in units of the predicted linear growth 
time, The dotted line represents the predicted linear growth of the instability. The 

transverse velocity is normalized to the initial perturbation. 

Fig. 7 Grey scale Images of the MHD Kelvin-Helmholtz instability. These images were taken from 
the high resolution (512 X 512) simulation at t = 1.4/T. Quantities shown (from left to right, 
top to bottom): gas density; magnetic pressure; gas pressure; magnetic field lines (contours 
of vector potential Az). 

Fig. 8 Gray scale images of gas pressure (upper left), magnetic pressure (upper right), V v (lower 
left), and (V X v) z (lower right) in the compressible Orszag-Tang vortex test. White indicates 
the regions with high (or positive) values and black indicates the regions with low (or negative) 
values. The calculation has been done in a periodic box with x = [0, 1] and y = [0, 1] using 
256 X 256 cells and a Courant constant of 0.6. The initial configuration has been set with p = 
25/367T, p = 5/127T, v = - sin(27ry)a; + sin(27rx)y, B = [- sin(27ry)a; + sin(47rx)y] /\/4tt, 
and 7 = 5/3 and the plot shown is at t = 0.48. 



